#####################
# Sleep Project - Pedro Bessone, Gautam Rao, Heather Schofield, Frank Schilbach, and Mattie Toma
# Purpose: Provides adjusted p-values for table 18 from the Appendix 
#          (Comparing Different Multiple-Hypothesis Testing Corrections for Table IV). 
# Last edited: 07 May 2021
#####################

rm(list = ls())

options(scipen = 99)
# install.packages("BiocManager")
require(multtest)
require(rio)
require(tidyverse)

# Load unadjusted P-Vals and divide them by family of correction 
df.p = import("Datasets/pvals_unadjusted.csv") %>% as_tibble()
colnames(df.p) = colnames(df.p) %>% str_remove('[0-9]')

df.family = df.p %>% 
  filter( var_name %in%  c("earnings", "wellbeing", "cogindex", "pref") ) %>% 
  mutate(order = seq(1,4)) %>% 
  relocate(order)

df.work = df.p %>% 
  filter( var_name %in%  c("prod", "typing", "output"))

df.wb = df.p %>% 
  filter( var_name %in%  c("physical", "psych"))

df.cog = df.p %>% 
  filter( var_name %in%  c("cogfunction", "attention"))

df.pref = df.p %>% 
  filter( var_name %in%  c("timepref", "social", "risk"))

# Load step-down FWER (Anderson 08) adjusted p-vals

df_step_down = import("Datasets/pval_adj_table_pool.csv") %>% as_tibble()

df_step_down_family = df_step_down[1:12,]

# Parameters
# corrections
corrections = c("BH", "BY", "hochberg", "holm", "TSBH")
alpha = 0.05 #(for rejection threshold when using the BKY06 method)

p.adjust.vec = function(p.unadjust, corrections, alpha){
  p.adj.mat = matrix(NA, 
                     ncol = length(corrections), 
                     nrow = length(p.unadjust))
  colnames(p.adj.mat) = corrections
  
  for(i in seq(1, length(corrections))){
    corr = corrections[i]
    
    if(corr == "TSBH"){
      ls.p.adj = multtest::mt.rawp2adjp(p.unadjust,
                                        proc = corr,
                                        alpha = 0.05)  
      p.adj.mat[,i] = ls.p.adj$adjp[,2]
    } else{
      ls.p.adj = p.adjust(p.unadjust, method = corr )
      p.adj.mat[,i] = ls.p.adj
    }
    
    
  }
  
  return(p.adj.mat)
} 



# Main table 


# Comparison: NS vs. Control

df.family.ns.cont = df.family %>% 
  select(order, var_name, p_val_ns) %>% 
  mutate(comparison = "NS vs. Control",
         order_2 = 1) %>% 
  rename(p_val_unadjust = p_val_ns) %>% 
  arrange(p_val_unadjust) 

p.unadjust = df.family.ns.cont$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.family.ns.cont = df.family.ns.cont %>% 
  bind_cols(p.adj.mat)

# Comparison: Nap vs. Control

df.family.nap.cont = df.family %>%
  select(order, var_name, p_val_nap) %>%
  mutate(comparison = "Nap vs. Control",
         order_2 = 2) %>% 
  rename(p_val_unadjust = p_val_nap) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.family.nap.cont$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.family.nap.cont = df.family.nap.cont %>% 
  bind_cols(p.adj.mat)

# Comparison: Nap vs. NS

df.family.nap.ns = df.family %>% 
  select(order, var_name, p_val_d_nap_ns) %>% 
  mutate(comparison = "Nap vs. NS",
         order_2 = 3) %>% 
  rename(p_val_unadjust = p_val_d_nap_ns) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.family.nap.ns$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.family.nap.ns = df.family.nap.ns %>% 
  bind_cols(p.adj.mat)

df.family.all = df.family.ns.cont %>% 
  bind_rows(df.family.nap.cont) %>% 
  bind_rows(df.family.nap.ns) %>% 
  relocate(var_name, comparison) %>% 
  rename(BKY06 = TSBH) %>% 
  arrange(order_2, order) %>% 
  mutate(ri_step_down = df_step_down_family$p_val_adj) %>% 
  select(!contains("order")) %>% 
  relocate(var_name, comparison, p_val_unadjust, ri_step_down, hochberg, holm, BH, BY, BKY06)

export(df.family.all, "Output/Appendix/Tables/TableA18/pval_adj_family.xlsx")


# Work Family


# Comparison: NS vs. Control

df.work.ns.cont = df.work %>% 
  select(var_name, p_val_ns) %>% 
  mutate(comparison = "NS vs. Control") %>% 
  rename(p_val_unadjust = p_val_ns) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.work.ns.cont$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.work.ns.cont = df.work.ns.cont %>% 
  bind_cols(p.adj.mat)

# Comparison: Nap vs. Control

df.work.nap.cont = df.work %>% 
  select(var_name, p_val_nap) %>% 
  mutate(comparison = "Nap vs. Control") %>% 
  rename(p_val_unadjust = p_val_nap) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.work.nap.cont$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.work.nap.cont = df.work.nap.cont %>% 
  bind_cols(p.adj.mat)

# Comparison: Nap vs. NS

df.work.nap.ns = df.work %>% 
  select(var_name, p_val_d_nap_ns) %>% 
  mutate(comparison = "Nap vs. NS") %>% 
  rename(p_val_unadjust = p_val_d_nap_ns)  %>% 
  arrange(p_val_unadjust)

p.unadjust = df.work.nap.ns$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.work.nap.ns = df.work.nap.ns %>% 
  bind_cols(p.adj.mat)


df.work.all = df.work.ns.cont %>% 
  bind_rows(df.work.nap.cont) %>% 
  bind_rows(df.work.nap.ns) %>% 
  relocate(var_name, comparison)

df.work.all = df.work.all %>% 
  rename(BKY06 = TSBH) %>% 
  arrange(comparison, var_name)

export(df.work.all, "Output/Appendix/Tables/TableA18/pval_adj_work.xlsx")


# Well-Being Family


# Comparison: NS vs. Control

df.wb.ns.cont = df.wb %>% 
  select(var_name, p_val_ns) %>% 
  mutate(comparison = "NS vs. Control") %>% 
  rename(p_val_unadjust = p_val_ns) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.wb.ns.cont$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.wb.ns.cont = df.wb.ns.cont %>% 
  bind_cols(p.adj.mat)

# Comparison: Nap vs. Control

df.wb.nap.cont = df.wb %>% 
  select(var_name, p_val_nap) %>% 
  mutate(comparison = "Nap vs. Control") %>% 
  rename(p_val_unadjust = p_val_nap) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.wb.nap.cont$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.wb.nap.cont = df.wb.nap.cont %>% 
  bind_cols(p.adj.mat)

# Comparison: Nap vs. NS

df.wb.nap.ns = df.wb %>% 
  select(var_name, p_val_d_nap_ns) %>% 
  mutate(comparison = "Nap vs. NS") %>% 
  rename(p_val_unadjust = p_val_d_nap_ns)  %>% 
  arrange(p_val_unadjust)

p.unadjust = df.wb.nap.ns$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.wb.nap.ns = df.wb.nap.ns %>% 
  bind_cols(p.adj.mat)


df.wb.all = df.wb.ns.cont %>% 
  bind_rows(df.wb.nap.cont) %>% 
  bind_rows(df.wb.nap.ns) %>% 
  relocate(var_name, comparison)

df.wb.all = df.wb.all %>% 
  rename(BKY06 = TSBH) %>% 
  arrange(comparison, var_name)

export(df.wb.all,  "Output/Appendix/Tables/TableA18/pval_adj_wb.xlsx")


# Cognition Family


# Comparison: NS vs. Control

df.cog.ns.cont = df.cog %>% 
  select(var_name, p_val_ns) %>% 
  mutate(comparison = "NS vs. Control") %>% 
  rename(p_val_unadjust = p_val_ns) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.cog.ns.cont$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.cog.ns.cont = df.cog.ns.cont %>% 
  bind_cols(p.adj.mat)

# Comparison: Nap vs. Control

df.cog.nap.cont = df.cog %>% 
  select(var_name, p_val_nap) %>% 
  mutate(comparison = "Nap vs. Control") %>% 
  rename(p_val_unadjust = p_val_nap) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.cog.nap.cont$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.cog.nap.cont = df.cog.nap.cont %>% 
  bind_cols(p.adj.mat)

# Comparison: Nap vs. NS

df.cog.nap.ns = df.cog %>% 
  select(var_name, p_val_d_nap_ns) %>% 
  mutate(comparison = "Nap vs. NS") %>% 
  rename(p_val_unadjust = p_val_d_nap_ns) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.cog.nap.ns$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.cog.nap.ns = df.cog.nap.ns %>% 
  bind_cols(p.adj.mat)


df.cog.all = df.cog.ns.cont %>% 
  bind_rows(df.cog.nap.cont) %>% 
  bind_rows(df.cog.nap.ns) %>% 
  relocate(var_name, comparison)

df.cog.all = df.cog.all %>% 
  rename(BKY06 = TSBH) %>% 
  arrange(comparison, var_name)

export(df.cog.all,  "Output/Appendix/Tables/TableA18/pval_adj_cog.xlsx")


# Preferences Family


# Comparison: NS vs. Control

df.pref.ns.cont = df.pref %>% 
  select(var_name, p_val_ns) %>% 
  mutate(comparison = "NS vs. Control") %>% 
  rename(p_val_unadjust = p_val_ns) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.pref.ns.cont$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.pref.ns.cont = df.pref.ns.cont %>% 
  bind_cols(p.adj.mat)

# Comparison: Nap vs. Control

df.pref.nap.cont = df.pref %>% 
  select(var_name, p_val_nap) %>% 
  mutate(comparison = "Nap vs. Control") %>% 
  rename(p_val_unadjust = p_val_nap) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.pref.nap.cont$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.pref.nap.cont = df.pref.nap.cont %>% 
  bind_cols(p.adj.mat)

# Comparison: Nap vs. NS

df.pref.nap.ns = df.pref %>% 
  select(var_name, p_val_d_nap_ns) %>% 
  mutate(comparison = "Nap vs. NS") %>% 
  rename(p_val_unadjust = p_val_d_nap_ns) %>% 
  arrange(p_val_unadjust)

p.unadjust = df.pref.nap.ns$p_val_unadjust

p.adj.mat = p.adjust.vec(p.unadjust, corrections = corrections, alpha = alpha) %>% as_tibble()

df.pref.nap.ns = df.pref.nap.ns %>% 
  bind_cols(p.adj.mat)

df.pref.all = df.pref.ns.cont %>% 
  bind_rows(df.pref.nap.cont) %>% 
  bind_rows(df.pref.nap.ns) %>% 
  relocate(var_name, comparison)

df.pref.all = df.pref.all %>% 
  rename(BKY06 = TSBH) %>% 
  arrange(comparison, var_name)

export(df.pref.all, "Output/Appendix/Tables/TableA18/pval_adj_pref.xlsx")
